Effects of fermions on the superfluid-insulator phase diagram of the Bose-Hubbard 

model 
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Building on the work of Fisher et al. (Phys. Rev. B 40, 546 (1989)), we develop the perturbation 
theory for the Bose-Hubbard model and apply it to calculate the effects of a degenerate gas of 
spin-polarized fermions interacting by contact interactions with the constituent bosons. For the 
single-band Bose-Hubbard model, we find that the net effect of the screening of the boson on-site 
interaction by the fermions is to suppress the Mott-insulating lobes in the Bose-Hubbard phase 
diagram. For the more general multi-band model, we find that, in addition to the fermion screening 
effects, the virtual excitations of the bosons to the higher Bloch bands, coupled with the contact 
interactions with the fermions, result in an effective increase (decrease) of the boson on-site repulsion 
(hopping parameter). If the higher-band renormalization of the boson parameters is dominant over 
the fermion screening of the interaction, the Mott insulating lobes in the Bose-Hubbard phase 
diagram are enhanced for either sign of the Bose-Fermi interactions, consistent with the recent 
experiments. 

PACS numbers: 67.60.Fp, 03.75.Mn, 03.75.Lm 



I. INTRODUCTION 

The superfluid to insulator quantum phase transition 
in a degenerate gas of bosons moving in a periodic po- 
tential is described by the Bose-Hubbard model intro- 
duced by Fisher et al. jl| almost two decades ago. In 
its simplest form, the model includes the hopping term, 
t, which describes the nearest neighbor tunneling ampli- 
tude of the constituent bosons on the lattice, and the 
on-site repulsion term, U, which approximates the short- 
range part of the Coulomb interaction if the bosons are 
charged (contact interaction if the bosons are neutral). 
In addition, the model includes the boson chemical po- 
tential, fi, which couples to the on-site charge density. 
The model can be directly implemented using spinless 
bosonic atoms moving in an artificially created periodic 
optical latticci^!^ 

It is remarkable that a theory developed based on 
the above simple premises can describe a true quantum 
phase transition with enough predictive power that can 
be tested experimentallyiii^ For large repulsive interac- 
tion U, boson charge fluctuations are suppressed and the 
system is in an insulating state. On the other hand, when 
the on-site repulsion is reduced, or, more appropriately, 
for large the system is in a superfluid state due to 
the Bose condensation of the mobile bosons. At some 
intervening value of jj, then, there should be a quantum 
phase transition separating the two phasesjii^ Recent ex- 
periments^i^iii^ using ultra cold bosonic atoms confined 
to an optical lattice, which mimics the model for periodic 
external potential in a custom setting, provided a first 
real demonstration of the superfluid to insulator transi- 
tion in an experimental system. The advantage of the 
atomic system lies in the ability to tune the parameters 
t, U and at will in a pristine, disorder-free environment. 



Even for periodic, disorder-free, external potentials, in 
a real solid state system additional fermions are always 
present and are invariably coupled to the constituent 
bosons. For example, in a granular superconductor,— 
where the Cooper pairs can be modeled as the bosons, 
there can be thermally generated quasiparticles.^" The 
question of additional fermions is also important in the 
context of the He3-He4 mixturesj^iii^ and the quark mat- 
ter, where two (color) quarks form a Cooper pair which 
interacts with the remaining unpaired quarks^ Remark- 
ably, such an additional degenerate gas of fermions can 
also be artiflcially introduced and coupled to the bosons 
in the ultra-cold atomic system.^** This raises an im- 
portant theoretical question as to what happens to the 
phases of the original Bose-Hubbard model^ when these 
additional fermions are present. The theoretical answer 
to this question can be experimentally tested in the so- 
called Bose-Fermi mixtures already realized in the optical 
lattice systems ji^ii^iii 

The effects of a degenerate gas of fermions 
on the Bose-Hubbard model have recently 
been investigated by various theoretical meth- 
ods >^>^i^>^>^>^>^'^i^>^i^>^i^>^^^i^>^i^>^ At 
first glance, the interaction between the bosons and the 
fermions seems to give rise to an effectively reduced 
repulsive interaction among the bosons compared to 
the bare (without the fermions) model. As a result, it 
may be expected that the superfluid phase coherence 
would increase in the Bose-Fermi mixture compared 
to the purely bosonic case. Recent theoretical and 
numerical works in Refs. 20ll3^ seem to agree with 
this preliminary assertion. Quite surprisingly, however, 
the opposite effect - the reduction of the superfluid 
coherence - was observed in the experiments via the 
measurement of the visibility of the quasimomentum 
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distributionJ^ii^iii The numerical work of Ref. [S^] 
argues that the decrease of the bosonic phase coherence 
may actually be due to finite temperature effects;^ 

The question of the intrinsic effect of the fermions on 
the Bose-Hubbard model is, however, far from theoreti- 
cally settled. The assertion of the reduction of the repul- 
sive interaction among the bosons is based on an hypoth- 
esis of static screening due to the fermionsr^-''^^ It has 
been recently argued in Ref. [36| that the dynamic part 
of the fermionic screening is equally important for the 
superfluid phase coherence of the bosons. By developing 
a Weiss-type, local, self-consistent mean field theory for 
this screening interaction, Ref. fSBj claims that an effect 
akin to the fermionic orthogonality catastrophe,"^^ aris- 
ing from the fermionic dynamic screening fluctuations, 
can suppress superfluidity. This way, the net intrinsic 
effect of the fermions may be in the same direction as in 
the experiments after all. 

In this work, we develop a rigorous perturbation the- 
ory to calculate the effects of an additional interaction 
potential (which, in the present context, is fermion- 
mediated) among the constituent bosons in the Bose- 
Hubbard model. To simplify the calculations, we take 
the system to be spatially homogenous, that is, we ne- 
glect the effects of the external confining potential in 
the optical set up. We also assume that the superfluid 
and the Mott insulating states are the only two pos- 
sible states, and neglect the possibility of other exotic 
states3>^>^>^>^>^>^>^ With these simplifying assump- 
tions, we take the largest interaction to be the Hub- 
bard U, and treat the additional fermion mediated in- 
teraction in perturbation theory. Using the Hubbard- 
Stratanovich transformation, we first rewrite the parti- 
tion function of the model in terms of a space- and time- 
dependent complex scalar field theory,^'^ whose coupling 
constants are given by the correlation functions of the 
original Hamiltonian modified by the fermion-mediated 
interaction. In mean field theory, the coupling constant 
of the quadratic term of the field theory, which is given 
by the boson on-site Green's function, provides the phase 
boundary between the Mott-insulating and the super- 
fluid phases.'* In the presence of the fermions, the boson 
Green's function must include the effects of the addi- 
tional fermion-mediated interaction. The calculation of 
the phase boundary is thus reduced to the perturbative 
evaluation of the boson Green's function in the presence 
of the additional space-time-dependent interaction. 

The above method still leaves us with a non-trivial 
problem, because the bare Bose-Hubbard model, on 
which we build the perturbation expansion of the Green's 
function in powers of the additional interaction, is not 
Gaussian (quadratic in the boson operators). As a re- 
sult, the standard machinery of bosonic perturbation 
expansiouf^S e.g.. Wick's theorem and the hnked cluster 
theorem which enable one to calculate the higher order 
corrections to the Green's function in terms of the in- 
tegrals over products of the bare Green's function, do 
not apply. Thus, one has to calculate the higher-order 



correlation functions non-perturbatively with respect to 
the Hubbard U. Fortunately, we need these correlation 
functions computed only with respect to the on-site part 
of the Bose-Hubbard model, which conserves the number 
of particles on every site. Because of this local number 
conservation, we are able to calculate these correlation 
functions exactly in the particle number basis {rii}. We 
also confirm that the apparent divergences, arising out 
of the summations over all the lattice sites and integrals 
over the imaginary time in the calculations of the cor- 
relation functions, are exactly canceled, and in the fi- 
nal result, one obtains non-zero perturbative corrections 
to the Green's function. Using this perturbatively cor- 
rected Green's function, we can calculate the effects of 
the fermions on the superfluid-insulator phase diagram. 

Using the methods outlined above, we find that, for 
the single-band Bose-Hubbard model fsections HIl IIIIl HVl 
ly]) . the fermions intrinsically shrink the area occupied 
by the Mott insulating lobes (Fig. [2]). The overall ef- 
fect is qualitatively in the same direction as in the effects 
of Ohmic dissipation in enhancing the superconducting 
phase coherence in Josephson junction array,^^!^''^ or in 
granular superconductors.^ This result is contrary to the 
orthogonality catastrophe argument of Ref. 1 36] , while it 
agrees with the numerical results of Ref. [32]. Experi- 
ments, however, have quite convincingly shown that the 
fermions expand the area occupied by the Mott insu- 
lating lobes. The earlier experiment a^^'^^ observed the 
loss of superfluid coherence for fixed attractive boson- 
fermion interactions, Ubf, which were larger in magni- 
tude than the boson on-site repulsion itself. Recently, 
this finding has also been confirmed for both attractive 
and repulsive interspecies interactions in a range of val- 
ues for 1?7bf1 both smaller and larger than [/ji^ From our 
rigorous perturbation theory, therefore, we conclude that 
the single-band Bose-Hubbard model is inadequate to ex- 
plain the experimentally observed loss of superfluidity of 
the bosons by adding a degenerate gas of fermions. 

Next, we treat the more general multi-band Bose- 
Hubbard model in the presence of the fermions (section 
IVI[) in the same analytical framework developed for the 
single-band model. Here, we first find that there is an ad- 
ditional effect on the bosonic system, due to the fermion- 
boson contact interactions, which is mediated by virtual 
transitions of the bosons to the higher Bloch bands of 
the multi-band model. This effect leads to an effective 
increase of the boson on-site repulsion, U, and a decrease 
of the hopping parameter, i, for either sign of the fermion- 
boson interactions. There is some numerical evidence of 
this effect (termed self-trapping) in Ref. 42] , for the case 
of attractive interspecies interactions only. We treat the 
two disparate effects - fermionic screening and the effects 
of the higher bands - within the same analytical frame- 
work. This theory provides a consistent explanation of 
the loss of bosonic superfluid coherence by introducing 
fermions, irrespective of the sign of the interspecies inter- 
actions, as seen in the recent cold atom experiments. As a 
bonus, the perturbation theory we develop for the Bose- 
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Hubbard model can be applied to calculate the effects 
of any additional interaction (not necessarily fermion- 
mediated) on the Bose-Hubbard phase diagram. For 
example, the effects of Ohmic dissipation^^ or the ef- 
fects of a second dilute gas of bosons on the superfluid- 
insulator phase diagram can also be evaluated by the 
methods described here. Some of our results were earlier 
presented in shorter forms in Refs. [30,3l|. We provide 
all the technical details of the theory in completeness in 
the current article, and mention that the results involving 
the renormalization of the bosonic hopping term due to 
the fcrmions in the presence of the multiband processes 
(section VI C), a mechanism of considerable quantitative 
importance in determining the final quantum phase di- 
agram, was not considered earlier by us and is thus a 
completely new result. 



II. SINGLE-BAND BOSE-HUBBARD MODEL 

The Hamiltonian for the single-band Bose-Hubbard 
model is given by. 



Hb — Hos + Hf , 



Hos^'^ ( ^fT-BiiflBi - 1) - P^BTlBi ) , 
i ^ 



(1) 
(2) 

(3) 



Here &J , bi are the spinless boson creation and annihila- 
tion operators on the site i and fiBi — h\hi is the boson 
density operator. [/ > and ^lb in the on-site part of the 
Hamiltonian, i?os, denote the on-site boson-boson inter- 
action and the bare chemical potential, respectively. The 
part of the Hamiltonian, iJt, that depends on the nearest 
neighbor pairs involves the nearest neighbor hop- 

ping matrix element, t^, for the bosons. The partition 
function for the model can be represented in terms of an 
imaginary-time path integral. 



where 



Zb = j S)6*S)bexp [-5b (6*, 6)] , 
SB{h*,h) dTh*drh + J drHE 



(4) 



(5) 



and P denotes the inverse temperature. 

By decoupling the boson hopping term using the 
Hubbard-Stratanovich transformation with a complex 
scalar field ipii'''): the partition function becomes, 

i 

X exp{-Sos[b*,h]+Sc[b*,h,i:*,tp,]). (6) 



Here, the symmetric matrix Wij has non-zero elements, 
tB, only for the nearest neighbors, (i, j). The on-site part 
of the action, 6*03, is given by, 

5os[6*,6,] = V / dTb*drb, (7) 

^ Jo 

+ J d.r'Y^ (^^'^^BiiflBi - 1) - f^BflBi^ ■ 

The coupled part of the action, Sc, is given by, 

[b* , 6. , ij* , = /''rfr^ (V* (r)6, (r) + (r)fe* (r)) . 



(8) 

To write the effective theory in terms of the scalar field 
V' only, we perform the cumulant expansion. 



zb = J ii^^.^r, e^v 



T 



(9) 



dTdr£[i/'(r , r), ip* (r, t)] 



where the Lagrangian 
/:[V'(r,r), V*(r,T)] = 



(10) 



Clip + C2 



dr 



Here we used the continuum limit for ipiir), i.e., ipiir) 
i]j{r, t). In this limit, the coefficient r is given by. 



ztb 



dT{TMr)b\m 



(11) 



where the brackets denote average with respect to the 
on-site part of the Hamiltonian i?os- Such imaginary- 
time-ordered averages can be conveniently calculated us- 
ing the path integral formalism with the on-site part of 
the action Sos- In mean field theory, r = gives the phase 
boundary between the insulator and the superfluid states. 
Thus, the problem of calculating the phase diagram of the 
model is reduced to the calculation of the one-particle, 
on-site, boson Green's function at zero Matsubara fre- 
quency G,(c^„ = 0), where G.(t-t') = -{TMT)b\{T')) . 
The phase boundary is then determined by r = 0: 



4 + G,(0)-0, 
ztb 



(12) 



where z is the coordination number of the lattice. The 
on-site Green's function can be easily calculated using 
the boson number basis {rii}, 



G.(T-r') 



e(r-r')e((5i;p)(no + l)e-("-"')^^'' 
+ Q{t' -T)Q{5Eh)noe^^-^'^^''A . (13) 
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Here Uq is the mean density of the bosons in the ground 
state at zero temperature, SEp = Utiq — ji and SE^ — 
jjL — U{nQ — 1) are the particle and hole excitation ener- 
gies, respectively. In the frequency domain, this function 



Here, the action S{b* , b, c' , c) is given by 



{no + l)Q[-ii + Uno] noe[fi-U{no~l)] 



iujn + ^i~Uno iuJn+fJ, — U{no — l) 

(14) 

Finally, using Eqs. and one arrives at the generic 
Bose-Hubbard phase diagram on the plane as 

shown by the solid curve in Fig. [2l 



III. SINGLE-BAND BOSE-HUBBARD MODEL 
WITH FERMIONS 

A. Partition function 

We consider a mixture of bosonic and spin-polarized 
fcrmionic atoms in an optical lattice. The full Hamilto- 
nian of the Bose-Fcrmi system is given hy H = Hb + 
Hp + Hbf, with Hp representing the fermionic part of 
the Hamiltonian and Hbf describing the inter-species 
interaction: 



Hp = -tp ^ (c\cj+H.c.\ - c|c 



HBF = UFB^nBi{clci - n°p^). 



(15) 
(16) 



, are the fermion creation and annihilation op- 
erators on site i, tp corresponds to the nearest neigh- 
bor hopping matrix element for the fermions, ^p is the 
fermion chemical potential, UpB describes the on-site 
boson-fermion interaction, and "nPp^ is the average den- 
sity of the fermions. In Eq. (|16p . the quantity mPp^ has 
been subtracted from the fermionic density, cjc^, to high- 
light the lowest order (in Upb) effect of the fermions on 
the constituent bosons, which is a trivial shift of the bo- 
son chemical potential: ^iB ^ l^B — UpB^fipi- Henceforth, 
this shift in the chemical potential is implicitly assumed 
in Hb- 

The partition function of the Bose-Fermi system is 
given as. 



j D6*D&DctS)cexp [-S{b\b, c^c)] 



(17) 



S{b*,b, c\ C) = ^ / dT{h* drh + c\drC,) + / (ItH. 

(18) 



To write an effective field theory analogous to that 
in Eq. (fTO]) we need to successively integrate out the 
fermions and the bosons as discussed below. 

B. Integrating out the fermions 



The first non-trivial effects due to the fermions ap- 
pear in the second order in Upb- By integrating out the 
fermions (note that the fermions appear only in quadratic 
order in iJ), the imaginary-time partition function be- 
comes. 



Z = Jj)b*me^p{-S,B[b*,b^]) (19) 

5eff [b* M]= dTj2 Kdrb^ + / dT (i^os + Ht) 

Jq j Jo 

-V" / dn dT2nBi(Ti)A'hj{Ti-T2)nBj{T2) 

(20) 



In the second order in U pb, the integral over the fermion 
degrees of freedom gives rise to an effective non-local 
density-density interaction for the bosons with the func- 
tion Mij{Ti — T2) defined as. 



My(ri -ra) 



FB 



{Anpi{Ti)AnFj{T2)) 



(21) 



In the frequency and momentum domain, the effective 
interaction Mq(Qn) is proportional to the fermion po- 
larization function. The exact form of the interaction 
Mq(ri„) depends on the dimensionality of the system. 
The effective fermion-mediated boson-boson interaction 
kernel in 3D is given by. 



^ FB 




1 - 



In 


'k- 




+ r 


1 


l-(^fc4 




2 

) 


In 






+ 1 






- 1 


^ 8fc 


" k , 




_fcH 




- 1 



(22) 



Here, Un = ^n/^Ep and k = q/2kF, with Ep and kp be- tively. A is the fermion mean level-spacing, A = l/vpV, 
ing the Fermi energy and the Fermi momentum, respec- 
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with Up the density of states at the Fermi level and V 
the volume of the miit cell, V = a'^. 



IV. BOSON GREEN'S FUNCTION IN THE 
PRESENCE OF FERMIONS 



C. Integrating out the bosons 

Using the Hubbard-Stratanovich transformation, we 
first decouple the boson hopping term, Ht, in Eq. (pO)) . 
We then integrate out the bosonic fields via cumulant 
expansion to find, 



(23) 



The calculation of the perturbative corrections to the 
bare boson Green's function, Eq. (fT5|) . is non-trivial be- 
cause the bare on-site Hamiltonian, Hqs, is not quadratic 
in the boson operators. Therefore, one cannot use the 
standard diagrammatic technique;^ because the Wick's 
theorem does not hold. To make progress, we write the 
average required for the corrected Green's function as. 
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Z = ZoJ S)^,DV'*exp(-5[V'„i/.*] 
where the action S[ipi, ip*] is given by. 



(24) 



{TMr)bmy - -y ^b*^bexp{-S'MT)b*{0) (28) 
= ^ J Sb*^bexp{-Sos + \SMT)b*{Q), 

where Z' is the partition function corresponding to the 
action S' in Eq. (|25p , which we have rewritten here using 
the definition 



In^cxp 



I dT^^bi{T)%};*{T) + H.c. \. = X! / ^'^i / dT2nBi{Ti)Mij{Ti-T2)nBj{T2 

Jo , / ,,, "'0 "'0 



The expectation value (...)' in Eq. (|24|) is taken with re- 
spect to the action Scs[b*,bi] with the boson hopping 
parameter ts — 0, i.e., with respect to the action 

S' ^ Sos-'Y] / d.Tl / dT2nBi{Tl)Mij{Tl~T2)nBj{T2). 

y Jo Jo 

(25) 

By expanding Slip, "0*] up to the fourth power of the field 
^p, and taking the continuum limit, we arrive at the action 
of an effective complex t/)^ field theory. 




with X — {r,r}. The coupling constants c'l, c'2, c' ,r' ,u' 
are given by the correlation functions of the bosonic fields 
with respect to the action S'. As before, the mean-field 
phase boundary between the superfluid and insulating 
phases can be obtained by setting the coefhcient r' to 
zero: 



(29) 

where A is used as a book-keeping parameter. Note that, 
from this definition, the linear order in A corresponds 
to the quadratic order in the boson-fermion coupling 
constant, Ufb, via Eq. (|2ip . Evaluating the perturba- 
tive corrections to the boson Green's function up to the 
quadratic order in Ufb, therefore, requires us to expand 
the second line of Eq. ((28)) up to the linear order in A. Ex- 
panding the numerator and the denominator of Eq. (|28p 
and collecting terms which are linear order in A, we get, 

{TMr)blm'={TMr)blm (30) 
+ X{TMr)bU0)Si)~X{Trb,{T)bl{0)){Si). 

Finally, putting A = f in Eq. pop . we derive, 
{TMr)blm'={TMr)blm (31) 

rP rP 

/ dTl dT2Mjl{Tl-T2)Kijl{T,Tl,T2). 
,, J-P J-0 



1 

zts 



rP 

/ drGUT)=0, 

J~3 



(27) 



where G-(t) — ~{Trbi{T)bl{0))' is the single-site boson 
Green's function, which, in the presence of the fcrmions, 
must now include the effects of the additional fermion- 
mediated density-density interaction. Thus, the problem 
is now reduced to the calculation of the on-site full boson 
Green's function by computing the corrections to Eq. (fT51) 
due to the fermion mediated interaction. As we show in 
the next section, this can be done perturbatively in Ufb- 



The higher-order correlation function Kiji(T,Ti,T2) can 
be conveniently calculated using second quantization rep- 
resentation (see Appendix [A]) . in which it is defined as 

Kiji (r, Ti , r2 ) = (r.r &j (t) 6j (0) rij (ti ) n; (r2 )) 

-{TMr)blm{Trn,{TMr2)). (32) 

Thus, we have reduced the problem of calculating the 
perturbative corrections to the boson Green's function 
due to the fermion-mediated, non-local interaction to cal- 
culating a higher order boson correlation function with 
respect to i/os given in Eq. 
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A. Higher order boson correlation functions 

Since the on-site part of the Bose Hubbard Hamil- 
tonian conserves the number of bosons, the correla- 
tion functions above can be calculated exactly using the 
particle- number eigenstates. From Eq. ([32)) . we see that 
there are two free spatial indices j, I and two free imagi- 
nary time indices ti, T2 in the definition of Kijiij, ri, T2). 
It is clear from Eq. (|3ip that the correction to the bare 
Green's function involves sums over the free spatial in- 
dices and integrals over the free imaginary times, both 
of which, in principle, can give diverging contributions 



(note that /3 ^ 00 as T — > 0). However, as we show 
in Appendix [XJ the second term on the right hand side 
of Eq. ([5^ subtracts these diverging terms exactly by 
constraining the spatial sums to i and limiting the ti , T2 
integrals to the imaginary time intervals between and r. 
In this way, we are able to get rid of the divergences in the 
perturbation theory, even though the usual linked cluster 
theoremSS does not apply to the interacting boson Hub- 
bard model. After subtracting the terms which would 
have produced divergent contributions in Eq. PT|) . the 
correlation function -/^^/(r, ri, T2) acquires the following 
form at T = 0: 



K,ji[T,Ti,T2) = e-'^-^Q[5Ep)Q{T){Q{T-Ti)Q{T-T2)Q{T2)Q{Ti) [{no + l)no{5,, + 5a) + {no + l)S^M (33) 
+ e(r2-r)e(r-Ti)e(T2)e(Ti)%7io(no + l) + e(ri-T)e(T-r2)e(T2)e(ri)^,,no(7io + l) 
+e(r-ri)e(-T2)e(Ti)5,,no(no + l) + e(r-T2)e(r2)e(-ri)(5,,no(no + l)} 

~e"''^^Q{5Eh)Q{-T) {Q{Ti-T)Q{T2)Q{-^i)nl6,, + Q{T2-T)Q{^T2)Q{Ti)nl5a +e{T2-r)Q{T-T^)Q{~T2)Q{^i)nl5a 
+ Q{Ti-T)Q{T-T2)<d{-T2)<d{-n)nl5,j 6 (ti - T )e (ts - r)e (-T2 )e (-Ti ) [n^ (5,, + 5,; ) - <5,j ,5,,7lo] } • 

I 



Here 5Ep and 5Eh are the particle and the hole excita- 
tion energies: 5E.p = Uuq — ^ and 5 Eh = /i — Uln^f — 1). 
It is important to note that the correlation function 
Kiji{T,Ti,T2) is irreducible and cannot be factored into 
the product of the bare Green's functions as would have 
been possible if Wick's theorem were applicable. 



B. Green's function in static approximation 

We now proceed to calculate the effects of the fermions 
on the boson Green's function in the static approximation 
(static limit of the fermion polarization function ri„ = 
and q ^ 0) neglecting the retardation effects of the 
induced interaction among the bosons. The expression 

for MqiQn) in Eq. ^ becomes Mq(0„) - 

We then substitute the corresponding expression for Mji , 

Mji{Ti —T2) — -j^SjiS{Ti — T2), into the correction term 
to the bare Green's function, which is the second term 
in Eq. (pij) . Because of the factor 6{ti — T2) in Mji{Ti — 
T2), we note that the correlation function Kiji{T,Ti,T2) 
in Eq. p3p reduces, in the static approximation, to 

K,jiiT,n,T2)~^ e(T)e(Ti)e(r2)e(T-Ti)e(T~r2) 

X [{Sij +Sii)nQ{no + l) + 6ijSii (no + l)] exp {-SEpT) 

+ e(-T)e(-ri)e(-T2)e(Ti - r)e(T2 - r) 

X [-(5,y + Suhl + StjSuno] exp (SEhr) . (34) 

Substituting this form of Kiji into the correction term 
to the bare Green's function, using the spatial Kronecker 
(5- function from the kernel Mji, and using the identity 



J2j I i^ij + Sii) — 2 J2j I ^ij 1 we find the correction in the 
static approximation to be given by, 

+ Oi-^) dTino(l-2no)e*-^"^^ . (35) 

Carrying out the ri integral and taking the Fourier trans- 
form to frequency space with respect to r, we finally de- 
rive the following expression for the full Green's function 
in the static approximation at zero frequency, 

U\no-^/U - (no - 1)/ 

_ / (no + l)(2no + l) _ no(2no - 1) \ 

2A(72 (no - ^i/U)^ - (no - 1))^ J " 

It is clear from this expression that near the degeneracy 
points {jj is an integer and tg — 0), where the gap in the 
Mott insulator states to the single particle excitations, 
5Ep/h , is small ~ Upg / A, the perturbation theory breaks 
down. This caveat would apply also to the calculation of 
the full Green's function given in the next subsection, 
and our results for the phase diagram would not be valid 
near the degeneracy points. 

In the static approximation, the same corrected 
Green's function as above could be obtained in an in- 
dependent, more obvious way. The method described 
above is useful, however, to calculate the corrections to 
the Green's function when the perturbing term to the 
Hamiltonian is explicitly dependent on space and imag- 
inary time, as in the case of the full fermion-mediated 
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FIG. 1: (color online) The dependence of the function R{y) 
on its argument. 



interaction. As we show below, the alternative, more 
transparent method to evaluate the Green's function in 
the static approximation produces the same answer as 
that found above, and this demonstrates the correct- 
ness of our perturbation formalism. In the alternative 
method, we could substitute the static, on-site form of 
Mij{Ti — T2) directly into the action, Eq. (PO)) . and calcu- 
late the Green's function exactly. It is easy to see that, 
in the static approximation, the mobile fermions simply 
renormalize fj, and U of the bare Bose Hubbard Hamilto- 
nian Hb: U - (7|b/^ and ^ ^ /z The 
exact Green's function, thus, can simply be obtained by 
substituting these renormalized parameters in Eq. (|14p . 
After expanding the result to the second order in Ufb, 
the resulting expression exactly matches that in Eq. ([55)1 . 



Full boson Green's function 



tion effect^ and the spatially non-local nature of the 
interaction kernel in Eq. (HI]). In order to take into ac- 
count these effects, we substitute the full expression for 
Mij{Ti — T2) into Eq. pTjl . Calling the second term on the 
right hand side of Eq. ((3T|) 6G[{t), we find the correction 
to the bare Green's function to be given by, 

SG',{t) = (e(T)(no + l)e-''^''M+e(-T)noe*^"(")) 

pT nT 

X dn dT2Mu{Tl-T2) (37) 

Jo Jo 

After performing the imaginary time integrals, we find 
that the required correction to the Green's function at 
zero frequency (see Eq. |2Z1)) is given by, 



13 



dT5G[{T) 



■/3 



E 

71— — OQ 



Mii{uJn)- 



dr siv? 



X [e(T)(no + 1)6-*^"" + e(-T)noe*^"n 
= 2 r ^M,,(a;„) 
'no + l 1 



no 



5Ej, ojI + SEI 5EhOjl+5El 



(38) 



Finally, using the form of A/q(f7„) from Eq. (1^^ . we 
obtain the following expression for the full boson Green's 
function at zero frequency, 



no 



1 



no 



U \nQ — ii/U fi/U — (no 



(39) 



AC/2 {{no-fi/Uy 
no 

{fi/U~{no-l)r^ 



R 



( " 




\AEf 


no 



u 

1e^ 



jj-{no-l) 



The static screening approximation for Afg(f2„) does 
not, however, take into account the important retarda- 



Here we introduced the dimensionless function R{y) 
given by, 



J 




~k 



In 



fc-if + 1 



In 



(40) 



where the momentum integral is defined in the First 
Brillouin zone with A — . The plot of the func- 
tion R{y) is shown in Fig. [1] As follows from Eq. ([39]) . 



the importance of the fermion renormalization effects 
is determined by the ratio of /x(~ U) and Ep. When 
the fermion density is small, i.e., fi/Ep ^ 1, the cor- 
rections to the Green's function are suppressed since 
R{y » 1) ^ 0. In the opposite limit, fi/Ep ^ 1, the 
function R(jj <^ 1) ~ 1, and thus, for a given value of 
UpB, the effects of the fermions on the bosons are more 



pronounced. 



SHIFT OF THE PHASE DIAGRAM WITHIN 
SINGLE-BAND MODEL 



As we discussed in section lllll using Eqs. (fT2l[T4l) which 
are applicable to the bare Bose-Hubbard model, one ar- 
rives at the bare Bose-Hubbard phase diagram on the 
(/i, ts) plane as shown by the solid curve in Fig. [51 Here 
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FIG. 2: (color online) Phase diagram of the Bose-Hubbard 
model with and without the fermions in 3D for the boson 
density no — 1. The phase diagram in 2D^ is qualitatively 
similar. Solid line describes the insulator-superfiuid phase 
boundary without the fermions. The dashed line corresponds 
to the same phase boundary when the effects of the fermions 
are taken into account. For comparison, we also show the 
dash-dot line which denotes the phase boundary when the 
effects of the fermions are only treated in the static approxi- 
mation, see text for details. The region near the degeneracy 
points (integer jj-), where the calculations of this paper do 
not apply, are implicitly excluded from this figure, see also 
discussion following Ea. (l36|l . Here we used, for illustrative 

purposes, x^— = 0.1 and -^j^=0.15. 



given by the dash-dot line in Fig. [2] One can see that, in 
the static approximation, the fermions markedly shrink 
the area of the Mott-insulating lobes in the phase dia- 
gram. Finally, using Eq. and Eq. ([55]) . we calculate 
the true phase boundary, within the single-band Bose- 
Hubbard model, as shown by the dashed line in Fig. [2l 
It is clear from this plot that the dynamic part of the 
fermion screening function indeed suppresses superfluid- 
ity, as argued in Ref. [S^. However, the net effect of the 
fermions, in a generic region away from the degeneracy 
points (where our calculations do not apply), is to still 
suppress the Mott-insulating lobes and enhance the area 
occupied by the superfiuid state. The sign of this overall 
effect is qualitatively the same as in the effects of Ohmic 
dissipation in enhancing the superconducting phase co- 
herence in Josephson junction array or in granular 
superconductors.— The single-band Bose-Hubbard model, 
therefore, proves to be inadequate in describing the loss 
of superfiuid coherence by adding fermions as seen in 
the experiments. In the next section, we will extend the 
analysis to the multi-band model in search of a consistent 
explanation of the experiments. 



VI. MULTI-BAND BOSE-HUBBARD MODEL 
WITH FERMIONS 

A. Effective three-band model for Bose-Fermi 
mixtures 



we have plotted the figure for 3D, while the figure in two 
dimensions'^" is qualitatively the same. Even though, in 
this figure, we have shown the phase boundary between 
the Mott insulating and the superfiuid phases for only 
the mean ground state boson density uq = 1, the results 
for other integer boson densities are qualitatively the 
same,^ and we have omitted them for simplicity. The 
Mott insulating states are characterized by r > 0, (V') — 
and survive in lobes extending from one degeneracy point 
to the next. The superfiuid state, on the other hand, is 
characterized by a non-zero value of the order parameter: 
r < 0, (-0) ^ 0. The transition between the Mott insu- 
lating and the superfiuid states is a continuous quantum 
phase transition which can be described by the 0^ field 
theory given in Eq. (fTO|) . 

In the static approximation for the screening, only the 
instantaneous part of the fermion polarization function is 
taken into account. The corresponding phase boundary is 



As we have seen in the last section, the effects of 
fermions on the phase diagram of the single-band Bose- 
Hubbard model is in a direction opposite to that seen in 
the experiments. It is then clear that the effects of the 
higher boson Bloch bands, which are so far neglected in 
the single-band model but can be significant in the ex- 
periments, should be taken into account. For the more 
general multi-band Bose-Hubbard model, there is an ad- 
ditional, higher-band, effect of the fermion contact inter- 
actions, which leads to an effective increase of the boson 
on-site repulsion and a decrease of the hopping param- 
eter. As we show below, when this higher-band effect 
dominates over the fermion screening of the bosonic inter- 
actions, it provides an explanation for the loss of bosonic 
superfiuid coherence by introducing fermions irrespective 
of the sign of the interspecies interactionsJ^ii^iii 

To elucidate the effects of the higher boson Bloch 
bands, let's start with the following second quantized 
Hamiltonian: 



H 



2mB 



M + Viat(r-) 



2m F 



- + Mat(r) 



$(r) 4 
l'(r) 



2 

9BF 



I d<ir$t(r)$t(^)$(r)$(r.) 

d'^r$^(a;)$(r)^'t(r)«'(r) 



(41) 
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Here V\at{r) denotes the lattice potential, V\iit{r) = 
^osin^(7r^) with a the lattice spacing. The cou- 
pling constants for the interactions are given by gsB ~ 
^7^^, Obf = ^m^" ' "^^^ei'e "^s is the mass of a 
bosonic atom, rrired is the boson-fermion reduced mass 
and QBE jBF ai'G the boson-boson and boson-fermion 
scattering lengths, respectively. <i>(r) and ^'(r) are the 
boson and the fermion field operators, respectively. 

We now expand the boson field operators in the Wan- 
nier basis, $(r) = Yl,i,a '^^ij - ri)bi^a, where Wa{r - ri) 
are the eigenstates of the single particle Hamiltonian 
=fB + Mat(r-) with Tb = /i. The indices 

a (and all other Greek indices used below) in 3D de- 
note a — {ax,ay,az) with a^.y^z — 1,2,3 labeling 3D 
vibrational modes. We consider below the dynamics of 
the bosons moving in the lowest few Bloch bands. For 
Vq ^ Efi, where Efi = tt"^ /2mBa? is the recoil en- 
ergy and Vq is the strength of the lattice potential, the 
Wannier wave functions Waix) can be locally approxi- 
mated by the wave functions of the ath excited state 
of a harmonic oscillator. Here we assume that lowest 
lying bands are well separated from each other which 
allows one to uniquely define optimally localized Wan- 
nier functions^i^, i.e Wa{r) are real symmetric or an- 
tisymmetric wave functions decaying exponentially with 
r. The fermion fields ^'(r) — J^i'^ai'i' ~ i'i)ci^ai where 
the fermion Wannier wavefunctions Ua{r — ri) are cho- 
sen using the mean-field Hamiltonian for the fermions, 
Tp + Vpir) with the effective mean- field potential for 
the fermions Vpir) — Vut{i') + 2^pB{r). Here, PBir) — 
n'o\wi^i{r)\'^ with pBif) and uq being the average boson 
density per site and average boson number per site, re- 
spectively. Thus, the shapes of these functions within 
a unit cell depend on the sign of the interspecies inter- 
actions. Restricting the fermions to a single band,— we 
get the following multi-band Hamiltonian for Bose-Fermi 
mixtures. 



H = Hi + H'i+Ht+ iJ, 



BF 



Hp 



Hi^'^eanBi,a + ^ ^^^-nBis{nBi,i-l) 

i,a i 

t(")[6L6,,„+H.c; 

<y>,Q 

HbF— ^ ^FB ["-^ct — {flBi.a) 



(42) 
(43) 

(44) 
(45) 
(46) 
(47) 



Here the boson single-particle and hopping energies are 

£a = (Wi^alH^lWi^a) - M and t^"'' = ~{Wi^a\TB + 

Vut(f)\wj,a) , respectively. The boson-boson and boson- 
fermion interactions are given by the following overlap 
integrals: Ubb = gBBiwiXiWi^ilw^Xi'^hi) and Upg = 



cjBF /'^{wi.a, Ui\wi^a', Ui).The fcrmion energy and hopping 
are, eo = (ui|Tp -I- VF{r)\ui) — fip and tp = —{ui\Tp + 
VF{r)\uj), respectively. The matrix elements M^^c^a are 
defined as M^^ctA = ^^{wi^p.;wi^i\wi^cy;wi^x)- In the low- 
est band, Ubb = -^iiii- The prime in the summation 
sign on the right-hand-side of Eq. indicates that 

Mini is excluded from the sum. 

Instead of treating the full, complex, multi-band 
Hamiltonian(^i^i^ we consider, for simplicity, a three- 
band model for the bosons and treat the fermions within 
a single band. We also keep only the largest band-mixing 
terms in the bosonic part of the Hamiltonian.— The ef- 
fective three-band model is justified for large inter-band 
energy separation n = V^^ErVq. In Ref. [31j|, we show 
that the band-mixing processes coupled with the fermion 
contact interactions, within an effective two band model 
for the bosons, lead to an increase of boson-boson repul- 
sion. Here, using an effective three-band model, we will 
show that the band-mixing processes additionally lead 
to a significant renormalization of the bosonic hopping 
parameter as well. Thus, the three-band effective model 
is the minimal model that captures all the significant 
multi-band effects. The band-mixing processes involving 
the higher bands (a — 4, 5, ...) have the same qualitative 
effects on the low-energy Hamiltonian. Therefore, in the 
rest of the paper we restrict our analysis to the lowest 
three bands a = 1,2,3. 



B. Renormalization of boson-boson interaction 

The renormalization of boson-boson interaction ap- 
pears due to the presence of the band-mixing terms 
given by H[ in Eq. ([1^ . In addition to the scattering 
of bosons within the lowest band given by the energy 
Ub b , there are processes involving scattering of bosons to 
higher Bloch bands. The particle(s) promoted to higher 
Bloch bands have much higher energy set by f2, and, 
thus, can stay in these excited states for a short time 
~ l/ri. Because of a large gap UBB,f^,t, higher 

Bloch bands can be integrated out yielding the renormal- 
ization of the parameters of the single-band Bose-Fermi 
model. Within the effective three-band model, the dom- 
inant band-mixing processes correspond to scattering of 
two bosons from the first band to two bosons in the sec- 
ond band described by the amplitude M2211 and scat- 
tering of two bosons from the first band to the first and 
the third bands given by Mami^ In order to calculate 
corrections to the boson Hamiltonian due to the band- 
mixing terms i?/, we use a perturbation series in l/fl. 
The lowest order corrections to Ubb appear in the sec- 



ond order and are proportional to Mljn/fi and M^m/il.. 
However, these renormalizations are present in a pure 
bosonic system as well and are independent of Up b ■ Al- 
though these are important corrections, we neglect them 
here since they are not modified by the presence of the 
fermions, which is the focus of the present paper. We 
will assume that these corrections are already taken into 
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FIG. 3: (color online) Virtual processes involving second (a) 
and third (b) Bloch bands leading to the fermion- dependent 
renormalization of the boson-boson interaction. Large (blue) 
and small (red) circles represent bosonic and fermionic atoms, 
respectively. 



account, yielding the renormalizcd Bose-Hubbard param- 
eters Ubb and t. 

Let us now concentrate on the renormahzations of 
the Bose-Hubbard parameters which are proportional to 
Ufb- The lowest non-trivial effects appear at the third- 
order in a perturbation expansion: 



AH, 



^ {no\Hl\s){s\HBFW)W\H'i\no) , , 



(48) 



in the lowest band ((?^s,i) = no), and the state |s) de- 
notes intermediate states where one or two bosons are 
excited to the higher vibrational states. The explicit eval- 
uation of the matrix elements yields the result, 



AHi = ^nFiUBiinBi-l)^ 



y^j'iFinBi{nBi-l)^y^^ 

i 11^1 



[E^+E^y 

Mh,,{U^/B-U¥B) 
El 



where Ex = + \y + \z — i) with \x,y,z being the 

quantum numbers labeling vibrational states along x, y, z 
directions, respectively. Within the effective three-band 
model, the correction to the boson-boson interaction be- 
comes 



AHi = ^ flFiflBiiflBi - 1) 



X 3 



^'^ 112 ^2 F^B^ ^Pb) 



3M 



FB 



(49) 
^U'f'b)' 



4172 

^nFiflBiiflB 



16f72 

rii 



i-,2 '^^'^1113x(^Fb'° ^Fb) 
' 4^(2 



Here the state I no) denotes the state with all the bosons 



Here the factors of 3 come from the symmetry of the 
overlap integrals, i.e., Mii2^2, = ^Ai2„2„ = ^112,2,- 
The magnitude of the overlap integrals can be evaluated 
within harmonic approximation. In general, we find that 
these overlap integrals quickly decrease with the increase 
of the band number. Therefore, the largest contributions 
to the renormalization of the boson parameters come 
from the lowest excited vibrational states. Using the 
values of the overlap integrals, we obtain the following 
correction to the boson-boson interaction. 



AH, 



^Ubb 
16172 



nFiflBiiflBi — l) 



U 



2A 
FB 



BF- 



^ BB 



y^^nFifiBiifiBi-l) 



-U 



FB 



9 

64 



-u. 



FB) 



^UIb 



32rj 



2 /L^ 



flFiflBiiflBi-lYiU] 



FB 



( 27 



V 1024 



32"^^ 



Pl3 



FB) 



(50) 



Here the dimensionless parameter pia is given by 



= 1 - 



FB 



U 



= 1- 



FB 



(wi^i;ui\wiX,Ui 



Using Eg. ([SO)) , the renormahzed boson-boson interaction 
Ucs becomes 



f/eff ^UbbII 



^UbbUbf 



8172 



-uf 



P12 



9 

64 



no 
2 



To understand the effects of the fermions on UcS, we 
first consider the case of attractive interspecies interac- 

(51) tions, Ubf < 0. The sign of the correction depends on 
the sign of pia- Since boson and fermions attract each 
other, the fermion wavefunction u(r) becomes peaked at 
the center of a unit cell. Therefore, the overlaps of 
with the boson Wannier functions in the second, Wi^2, 
and third, w,;.3, Bloch bands are smaller than its overlap 
with Wi^i. Thus, for negative Ubf, P12 and P13 are pos- 

(52) itive and 0(1). Consequently, the presence of fermions 
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FIG. 4: (Color online) Virtual processes leading to the renor- 
malization of the bosonic hopping. Here the sequence of 
events corresponds to the amplitude A^^^ (see text). Large 
(blue) and small (red) circles represent bosonic and fermionic 
atoms, respectively. 



leads to the increase of the boson-boson repulsion. 

In the opposite hmit, Ubf > 0, the fermions and 
bosons repel each other within the unit cell. Assum- 
ing that the number of the bosons per site is larger than 
that of the fermions, the fermion density is likely to be 
suppressed at the center of a unit cell, resulting in the 
numerator in the second term in Eq. (|5ip exceeding the 
denominator, p < 0. Since the sign of the renormaliza- 
tion to UcS in Eq. (j52p is determined by sgii(pUFB)i the 
renormalization is again repulsive, and it boosts Ubb as 
in the case of Ub f < 



C. Renormalization of boson hopping parameter 

In addition to the renormalization of Ubb, the pres- 
ence of the fermions leads to a renormalization of the 
boson hopping energy as well. The higher-band renor- 
malizations of t appear in the fourth order of the pertur- 
bation theory, and, thus, in principle, are smaller in l/fl. 
However, the boson tunneling rate between the nearest 
neighbor sites is much larger for the second and third 
bands compared to the first one, 

^(1) ^ ^(2,3) _ Therefore, 
the higher-band corrections to the hopping energy can 
also be significant. 

Let us consider the tunneling process of a boson 
to a nearest neighbor site through the higher Bloch 
bands. Again, we find that there are significant correc- 
tions to hopping which are independent of the boson- 
fermion interactions. These virtual processes lead to 
the renormalization of the hopping parameter within 
the bare Bose-Hubbard model. As before, we ignore 
these processes below. We now consider the lowest-order, 
fermion-dependent, tunneling processes through the vir- 
tual states. The lowest-order corrections to the tunneling 



Hamiltonian are 
AHt^ ^ A{n,m){\ni,mj){n-li,m + lj\ +II.C.} 

(53) 



Here \ni,mj) is the state of the Mott insulator with n 
and m bosons on sites i and j, respectively. The state 
\n—li, m+lj) = \n— 1)^ \m + l)j where the nearest- 
neighbor sites i and j have one less and one more boson 
with respect to the original state \ni,mj), and the oc- 
cupation of all other sites remains the same. One can 
now calculate the lowest-order fermion-dependent ampli- 
tude A(n, m) for the tunneling process . Assuming that 
the initial state and final states are |i) — \ni,mj) and 
I/) = |n— li, m-|-lj), respectively, the amplitude A(n, m) 
is given by 

A(n,m)=A(i) + (54) 
mHi\s){s\H,\s'){s'\HBF\s'){s'\Hl\f} 



(55) 

,(2) „ V (»l HI \s){s\ Hbf \s){s\ Ht \s'){s'\ H[ \f) 

{E-E,){E-E,Y ■ ^ ^ 



Here \s) denotes the intermediate states with one or two 
excited bosons. The dominant contribution to the am- 
plitude A comes from the virtual processes involving the 
third Bloch band, see Fig ID Virtual processes involv- 
ing the second band are proportional to (i^^^)^ since 
they require to transfer two bosons from site i to site j. 
Therefore, their contribution to A{n, m) is much smaller 
compared to the one from the third band, and can be 
neglected. The resulting expression for the amplitude 
A(n, to) becomes 



A{n, to) : 



r!3 



7-r(l,l)_7-r(3,3) 
'^FB '-^FB 



(n — l)my^n{m + l) 



(57) 



Note that the amplitude for tunneling through virtual 
state is zero if either n = 1 or m = since the band- 
mixing processes require to have at least two bosons on 
the same site. In the Mott-insulating state the aver- 
age number of bosons per site is np. Hence, the domi- 
nant contribution to the hopping Hamiltonian (|53p comes 
when m — n = uq. (In the Mott-insulating state, one 
can neglect the fluctuations of boson occupation per site 
at the mean field level.) Using Eqs. (|53p and ([57| and 
evaluating the overlap integrals in the harmonic approx- 
imation, we find the renormalized boson hopping energy 
in the Mott-insulating state to be 



te«— t 



t^^^ ^BS ('^0 - l)"-o V "0 [no - 



1),, F 



(58) 
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with the function pi3 defined in Eq. 
matrix clement t*-^-* is given by, 



5T|) . The tunnehng 



5.0 r 



i(3) _ 

~ 4 



Ma 3 



A/a 2 



(59) 



where Afyi(r, x) is the characteristic Mathieu function.^ 
In Eq. ([58|) . the hopping energy for the lowest band 
i'^^^ includes all the corrections due to the fermion- 
independent virtual processes (fermion-independent 
higher-band correction to the hopping energy goes as 



t^'^^C/^^/ri^). Therefore, the self-consistent per- 



turbation theory for the tunneling is well-defined, and 
the second term in the brackets in Eq. ([58]) is smaller 
than one. However, given that UggllpB 1 and 
^(3)^{(i) ^j^g renormalization of the boson timncling 

is still a significant effect. 

We now discuss the sign of the correction to the boson 
hopping parameter. As mentioned before, in the case 
of attractive interaction between the fermions and the 
bosons {UpB < 0), the function pi^ > 0. Therefore, 
the higher-band virtual transitions lead to a suppression 
of the boson hopping parameter. This can be also un- 
derstood using the arguments in Ref. the presence 
of fermions leads to the squeezing of the bosonic Wan- 
nier wavefunction, and, thus, should reduce the bosonic 
hopping energy. In the case of repulsive inter-species in- 
teraction {UpB > 0), the bosons and the fermions like 
to maximize the distance between each other. Assum- 
ing, as a result, that the fermion density is substantially 
suppressed at the center of the unit cell, the function pi3 
becomes negative. Thus, the presence of the fermions 
leads to the suppression of the bosonic hopping energy 
in this case as well. 



D. Shift of the phase diagram 

We now discuss the effects of the fermions on the phase 
boundary between the Mott insulator and the superfluid 
states. Within the effective three-band model, we find 
that the parameters of the bare Bose-Hubbard model are 
renormalized as Ubb Ucs and i tcs according to 
Eqs. ((52)) and ([58]) . The calculation of the phase bound- 
ary can be done in a similar manner as in Sec|Vl In- 
tegrating out the fermions, and neglecting terms of the 
order of leads to the effective imaginary-time 



action analogous to Eq. ([20 



b* dr hi -f ^^flBi {flBi - 1 ) - Mbi 



5eff(6*,6)=V / dT 

Jo 

13 . 

dr^teff(6lfoj+6]&») (60) 

<ii> 

dri I dT2nBi{Ti)Mij{Ti-T2)nBj{T2). 

AA Jo Jo 
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FIG. 5: (Color online) Phase boundary of the boson Hubbard 
model in 3D with and without the fermions. Here we choose 
the boson density no = 5 for illustration purposes. Solid line 
describes the insulator-superfluid phase boundary without the 
fermions. The dashed line corresponds to the phase boundary 
in the presence of the fermions. Here we used the parameters 
specified below Eg. fTO)) and Ubf = IJbb/'2. 



Here [/oft = Ubb + SUbb, tcs ^ i + St, fl ^ ^ - UpBn^^, 
with n}p^ the average density of the fermions. The change 
in boson-boson interaction SU b b and bosonic hopping St 
are given by 



SU 



BB 



SU'^gllBF 









Pl2 + 







St = t 



(3) U]jB{no-^)nn^yno{no + l) 



4fi3 



UpBPisni 



(61) 



(62) 



The last term in Eq. (pD)) describes the screening of the 
bosonic repulsive interactions by the fermions, same as 
in the single-band model, leading to the suppression of 
the Mott insulating phase. 

The effect of the two competing contributions - 
fermion screening and the effects of the higher bands - 
on the phase diagram can be calculated analytically using 
the framework developed in the preceding sections. We 
again need to calculate the boson on-site Green's function 
for the action in Eq. ((SD|) at zero frequency. 



(63) 



"-0 + 1 

SE. 



p 



SUBBno 



TP 



SEp ASEp 



SEp 
AEf 



no 
SEh 



1 + 



SUBBjno - 1) 
SEh 



—JIM-R 

ASE,^ 



SEh 
4:Ef 



Here SEp and SEh are the new, renormalized, particle 
and hole excitation energies: SEp = Ubb^q — jl and 
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SEh = {1 — UBB{no — I), where no is the number of 
bosons per site minimizing the ground state energy. 

As before, the mean field superfluid-insulator phase 
boundary can be obtained by solving the new equation, 
analogous to Eq. (P7)) . 



1 



+ / drGr(T) - 



(64) 



The shift of the superfluid-insulator phase boundary can 
be obtained by looking at the change of the critical hop- 
ping tc where the transition occurs with and without 
fermions, 5tc = tdUsF) — tdUsF = 0). To the linear 
order in UpBi the correction to the phase boundary is 
given by. 



St 



(1) 



St 



su 



BB 



Ubb 



Ubb zUbb 



(65) 



[1 + 2{Ii/Ubb? - 2ifl/UBB){na - 1) - no]no 



[1 + WUbb 



Notice that, as discussed before, the product of Ubf 
and p remains negative, irrespective of the sign of Ubf- 
Therefore, St < and SUbb > 0. For /I/Ubb, we use the 
value near the tip of the Mott lobes in the bare model 
(without the fermions): ji/UBBltip — \/no{rio + 1) — 1. 
By substituting this value into Eq. ([65]) one finds that the 

shift St'^^ is given by 



St 



(1) 



Ubb 



St 



SU, 



BB 



tip 



Ubb zUbb 



(l + 2no - 2A/no(l + no)) . 

(66) 



One can see that the above is always positive, indicating 
that the Mott-insulating phase expands at the tip of the 
lobes. 

We now consider the effect of the fermion-mediated 
screening, which manifests only in the second order in 
Ufb- The shift of the critical hopping, Stc , is given by, 

St^^ 



U'fb 



Ubb 



zAU 



BB 



1 



Ubb 



(67) 



X {{no + 1) 



U. 



(no-1) R 



BB 



Ubb 
4Ef 



no- 



^BB 



+no no--=- 



Ubb 



R 



'Ubb 
AEf 



Ubb 



-(no-1) 



which indicates that St^c^ < 0, and, thus, the fermion- 
mediated screening leads to the suppression of the Mott 



f 2) 

insulating phase. At the tip of Mott lobes, Stc becomes 
Stf^ U'fb (l + 2no - 2v/no(l + no)) 



Ubb 



xiR 



zAUbb 
UBB[\/na{l+no)-no] 
^Ef 



(68) 
(69) 



, „/c/BB[l+no-v'no(l + no)M 

V )}' 

We now compare the two effects, multi-band and 
fermion-mediated screening, for no S> 1. At the tip of 
the Mott lobes, the two corrections to the phase bound- 
ary come with a ratio. 



St 



(1) 



St 



(2) 



-AzTLoSt + su. 



BB 



FB 

A 



Aznt 



!!!! £ 



n^\UBF\Ril§f) 



(70) 
(71) 



For the purposes of illustration, we choose some typical 
experimental parameters: no = 5, z = 6, n'p = 0.75, 
Er = 2.6 • 10-30J, Vo/Er = 9, n/ER = 6, Ubb/Er = 
0.2, t'^^'^/ER = 0.9, A/Er = 47 and Ep/Er = 0.8 • lO^^. 
Assuming Ubf < 0, one can estimate pi^ « 0.6. For 
these parameters the function Rij§f:) ~ 0.02, and the 

ratio \Sti^'>/Stl?^\ becomes \Sti^'> / St'^?^ \ - Ubb/\Ubf\-10^ 
indicating that for realistic parameters the higher-band 
effect is dominant. Thus, we conclude that the super- 
fiuid state is suppressed for either sign of the interspecies 
interaction, as shown in Fig. [51 



VII. SUMMARY AND CONCLUSION 

In this paper, we have addressed the quantum phase 
diagram of the Bose-Hubbard model in the presence of a 
degenerate gas of spin-polarized fermions with the bosons 
and the fermions interacting via contact interactions. We 
have addressed this question by developing a framework 
for carrying out perturbation theory in the Bose-Hubbard 
model. We first conclude that, for the single-band Bose- 
Hubbard model, the degenerate gas of fermions intrin- 
sically shrinks the area occupied by the Mott insulating 
lobes (Fig. [2]) . For the multi-band Bose-Hubbard model, 
which is more general and experimentally relevant, we 
show that the virtual transitions of the bosons to the 
higher Bloch bands, coupled with the contact interactions 
with the fermions, generate a new renormalization of the 
interaction and the hopping parameter of the bosons. For 
either sign of the coupling between the fermions and the 
bosons, this renormalization enhances the boson on-site 
repulsion and decreases the hopping parameter, favoring 
the Mott insulating phase. If this effect is dominant over 
the fermion mediated screening effect, the superfluid co- 
herence of the Bose-Hubbard system will be suppressed 
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by the fermions, irrespective of the sign of the interspecies 
interaction, as has been observed in recent experiments. 
-'^^'^^1^^ Thus we conclude that the reconciliation of the 
experimental observations with the theory of the Bose- 
Hubbard model lies in the higher-band effects, which is 
also supported by the numerical calculations in Ref. [13] . 
Note that, using our general analytical framework, we 
have been able to explain the loss of superfluid coher- 
ence for either sign of the Bose-Fermi interaction, while 
the calculations of Ref. [i^l apply only to the case when 
the interspecies interaction is attractive. 
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We have not discussed in this paper the fluctuation 
effects on the phase diagram , ^^'^'^i^'* and the fate of the 
universality class of the superfluid to insulator quantum 
phase transition in the presence of the fermionsi^ Based 
on tree level arguments, Ref. [5^ argues that the critical 
properties of the transition at the tip of the Mott insulat- 
ing lobes, where the dynamic critical exponent z — 1 (the 
coupling constant ci = in Eq. (|10p ). are strongly mod- 
ified by the fermions. On the other hand, for a generic 
point on the boundary of the lobes, where z — 2, the 
critical fixed points are stable with respect to the fermion 
mediated interactions. The actual fate of the z — I fixed 
point, however, is an open question, and to answer this 
one has to go beyond the tree level of the renormaliza- 
tion group analysis. It is also important to note that the 
fermions mediate a time- and space-dependent density- 
density interaction among the bosons. Such spatially 
non-local interactions may give rise to spatially ordered 
states, which may or may not coexist with superfluidity. 
If they do coexist with superfluidity, there may an inter- 
esting possibility of inducing bosonic supersolids in some 
parameter regime by tuning the boson-fermion interac- 
tion. 



In conclusion, we develop a perturbation theory for the 
Bose-Hubbard model, which is different from the stan- 
dard theory for bosonic systems,''^ but is appropriate 
for the Bose-Hubbard model with a quartic (Hubbard- 
U) term. Using this theory, and by including the effects 
of the higher boson Bloch bands in a multi-band formula- 
tion, we explain the observed, unexpected, expansion of 
the Mott insulating lobes of the Bose-Hubbard model by 
introducing fermions. Our three new important findings 
are: (1) There are independent physical contributions 
to the bosonic quantum phase diagram in the presence 
of fermions, namely, the modification of bosonic interac- 
tion due to screening and the interband virtual transition 
processes in the multiband case; (2) the multiband pro- 
cesses are independent of the sign of the interaction, and 
if dominant, would lead to suppression of superfluidity, 
as observed experimentally, independent of the sign of 
the fermion-boson interaction; (3) the renormalization of 
the bosonic hopping term by the multiband transitions 
appears to be the most important process quantitatively. 



APPENDIX A: CALCULATION OF THE 
CORRELATION FUNCTION K,ji{t,ti,t2) 

In this Appendix we illustrate how we calculate the 
correlation function Kiji{T,Ti,T2), and show that the 
terms that can give divergent contributions to Eq. ([51]) 
exactly cancel out. To illustrate the cancelation, we 
consider the correlation function Kiji{T,Ti,T2) defined 
in Eq. Given that the on-site part of the boson 

Hubbard Hamiltonian conserves the number of bosons, 
the correlation function Kiji{T,Ti,T2) can be calculated 
in the second quantized representation by inserting the 
representation of unity as a sum over particle-number 
eigenstates. Doing this, we find that the second term in 
Eq. ( [321) zero temperatures is given by 

{TMr)blm{Trn,iT,)ni{T2)) ^ (Al) 
= e-*^-"n2(,io + l)e(r)e((5£;p)+e*^''^7i3e(~r)e((5£;,,), 

where no, 6Ep and SEh were defined in the text. To cal- 
culate the first term in Eq. ([5^ . we note it is important 
to distinguish between the cases where the indices j and 
I are the same or different from the index i, because they 
give rise to different matrix elements. To carefully take 
this into account, we separate the sums over i, j as j = 

E^J [i^-6^3){l~Sa) + il- % )Sa + (1 - 6a)S,j +S,,6a]. ' Af- 
ter taking care of the imaginary time orderings between 
the time indices 0, ri, r2, t and calculating the matrix el- 
ements, we find that the first term of Eq. is given 
by, 

{TMr)bl{0)n,{n)ni{T2)) = (A2) 

i=1..12 i=1..12 

Here, the functions li are given by, 

h = e(T2-Ti)e(Ti-T)e(T2)e(ri)n2(„o+i) (as) 

/2 = e(Ti-T2)e(T2-T)e(T2)e(ri)n2(„o + l) 
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/3=e(T2-T)e(r-ri)e(T2)e(Ti)x 

X [nl{no + l) + 6ijno{na + l)] 
/4 = e(ri-r)e(r-T2)e(r2)e(Ti)x 

X [nl{no + l) + S^ino{no + l)] 
h = e(r-ri)e(T-r2)e(ri - T2)e(T2)e(ri) x 

X [nl (riQ + 1 ) + (no + l)no (% + (^i; ) + (»^o + 1 ^jz] 
/6 - e(r-ri)e(T-r2)e(r2 - Ti)e(T2)e(ri) x 

X [no (no + 1 ) + (no + l)no (5^ + (5^; ) + (no + l)(5y (5^;] 
h = e(r2-r)e(r-ri)e(r2)e(-Ti)ng(no + l) 
4 = e(ri-r)e(T-r2)e(-T2)e(Ti)n2(no + l) 
/g = e(r-ri)e(-r2)e(Ti) [no(no + l) + 5yno(no + l)] 
/lo = e(T-T2)e(T2)e(-ri) [no(no + l) + 5ji(no + l)no] 
/ii - e(-T2)e(Ti -T2)e(r)n2(no + l) 
/i2 = e(-T2)e(T2 -Ti)e(r)n2(no + l), 

and the functions Ji, which are of similar form but cal- 
culated for the other 12 time orderings corresponding to 
(?(— t), arc omitted here for simplicity. 

It is straightforward to check that the terms in 
Eq. (jA3p which are independent of the site indices i,j 
can be summed to yield nQ(no + 1)- Thus, by comparing 
Eqs. (|A1|) and (jA2p one can see that these terms, which 



would have given divergent contributions to the Green's 
function (see Eq. ((3T|) ) cancel out leading to the r > 
part of Eq. (|33l In a similar way, for r < 0, the terms 
in Ji which are independent of the site indices sum up 
to ng and exactly cancel with the corresponding term in 
Eq. (|X2|) . leading to the r < part of Eq. dSH]). Fur- 
thermore, since in the static approximation the effective 
interaction kernel Mij{Ti — T2) itself imposes ri = T2, only 
the time orders e(T)e(ri)e(r2)e(T - Ti)e(r - T2) and 
e(-T)e(-Ti)e(-T2)e(Ti - r)e(T2 - r) from Eq. ^ 
contribute to the correlation function and one ends up 
with Eq. (|34p . However, for these time orders, the ti,T2 
integrals are constrained to be in the interval between 
and r and do not diverge. In a similar way, all the 
integrals in the calculation of the Green's function for 
the full interaction kernel are restricted to finite inter- 
vals as well, leading to finite non-zero contributions. For 
the calculation of the boson Green's function with the 
full interaction kernel Mq(f2„) given in Eq. ([22|) . we can 
use the fact the it is zero for either q or Qn equal to 
zero. As a result the terms of Eq. (|33l) that contribute 
to Eq. ((3T|) have i = j — I and the same time orders as 
given above. This way all the integrals in the calculation 
of the Green's function are restricted to finite intervals, 
canceling the potential divergences. 
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